Comparison of Actinobacteria communities from human‐impacted and pristine karst caves

Abstract Actinobacteria are important cave inhabitants, but knowledge of how anthropization and anthropization‐related visual marks affect this community on cave walls is lacking. We compared Actinobacteria communities among four French limestone caves (Mouflon, Reille, Rouffignac, and Lascaux) ranging from pristine to anthropized, and within Lascaux Cave between marked (wall visual marks) and unmarked areas in different rooms (Sas‐1, Passage, Apse, and Diaclase). In addition to the 16S rRNA gene marker, 441 bp fragments of the hsp65 gene were used and an hsp65‐related taxonomic database was constructed for the identification of Actinobacteria to the species level by Illumina‐MiSeq analysis. The hsp65 marker revealed higher resolution for species and higher richness (99% operational taxonomic units cutoff) versus the 16S rRNA gene; however, more taxa were identified at higher taxonomic ranks. Actinobacteria communities varied between Mouflon and Reille caves (both pristine), and Rouffignac and Lascaux (both anthropized). Rouffignac displayed high diversity of Nocardia, suggesting human inputs, and Lascaux exhibited high Mycobacterium relative abundance, whereas Gaiellales were typical in pristine caves and the Diaclase (least affected area of Lascaux Cave). Within Lascaux, Pseudonocardiaceae dominated on unmarked walls and Streptomycetaceae (especially Streptomyces mirabilis) on marked walls, indicating a possible role in mark formation. A new taxonomic database was developed. Although not all Actinobacteria species were represented, the use of the hsp65 marker enabled species‐level variations of the Actinobacteria community to be documented based on the extent of anthropogenic pressure. This approach proved effective when comparing different limestone caves or specific conditions within one cave.


| INTRODUCTION
Limestone caves are isolated habitats that are extremely limited in organic carbon . They are associated with the development of microbial strategies that enable adaptation to oligotrophic and high calcium conditions. In addition, caves are considered to be extreme environments, as they are oligotrophic, due to the complete darkness, low and constant temperatures, and high humidity (Barton & Jurado, 2007;De Mandal et al., 2017;Ortiz et al., 2014). Some of these caves contain Paleolithic artwork and, consequently, have been frequently visited by tourists for a long time (Alonso et al., 2019;Bontemps et al., 2021;Schabereiter-Gurtner et al., 2002). Thus, human-mediated dispersion of endogenous as well as exogenous nutrients and microorganisms inside the caves, together with human metabolism may induce changes in internal microclimatic and microbiological conditions (Hoyos et al., 1998;Mulec, 2014). Microorganisms can be somewhat resilient to changes Tomczyk-Żak & Zielenkiewicz, 2016). However, community shifts may also occur, which deserves attention if those shifts favor microorganisms with biodeterioration ability (F. Bastian et al., 2010;Bontemps et al., 2021;Fernandez-Cortes et al., 2011;Sánchez-Moral et al., 1999) or increase the proportion of species harboring strains with pathogenic potential (Bercea et al., 2019;Neral et al., 2015;Rajput et al., 2012). In such cases, these shifts could represent a threat to cave conservation and the health of visitors.
However, some species, such as Pseudonocardia hispaniensis, have been detected in the moonmilk of a cave associated with urban contamination , and its presence may be due to its ability to degrade exogenous organic compounds (Lavoie et al., 2017;Porca et al., 2012). This is also true for the Streptomycetaceae family, which has also been identified among pristine cave inhabitants (Maciejewska et al., 2017). Furthermore, Pseudonocardiaceae and Streptomycetaceae seem to be prominent in pigment-forming communities on the walls of certain show caves (Cuezva et al., 2012;Porca et al., 2012).
In contrast, Euzebyales (Gonzalez-Pimentel et al., 2018), Nocardiaceae (De Mandal et al., 2017), and Gaiellales (Zhu et al., 2019) represent oligotrophic taxa that dominate in pristine caves, most likely because of their ability to degrade complex compounds present in the cave, or to fix CO 2 .
Identification of Actinobacteria taxa in these environments has been based on culture-dependent approaches and 16S rRNA gene sequencing (De Mandal et al., 2017;Yasir, 2018). Culture-dependent approaches limit community coverage and should, therefore, only be used as a supportive determinant to study their specialized metabolisms (Hamedi et al., 2019;Long et al., 2019). The 16S rRNA gene is the most frequently used marker for bacterial systematics; however, it is unable to resolve individual sequences to the species level (Fox et al., 1992;Hirsch et al., 2010). In addition, identical 16S rRNA gene sequences can be found in different taxa or different 16S rRNA gene sequences within a single species (Větrovský & Baldrian, 2013). However, other (protein-coding) genes with high sequence polymorphisms used for multilocus sequence analysis (MLSA; i.e., rpoB, gyrB, and hsp65) may help distinguish closely related species from environmental DNA. Thus, they may present greater discriminative power at the species level compared with the 16S rRNA gene (Gao & Gupta, 2012;Takeda et al., 2010;Vos et al., 2012). Thus, using a metabarcoding approach, Aigle et al. (2021) utilized another gene (tpm) to evaluate the biodiversity of the bacterial community associated with the degradation of pollutants in anthropized environments at a species level.
In the case of Actinobacteria, hsp65 encodes a molecular chaperone (a 65-kDa heat-shock protein), which, when amplified by specific primers (Telenti et al., 1993), has been used as a molecular marker to identify a large number of isolates (Rodríguez-Nava et al., 2006. A single copy of this gene is present in Actinobacteria genomes, unlike the 16S rRNA gene marker (Colaco & MacDougall, 2014). Additionally, the length of the amplified hsp65 fragment (441 bp) is sufficient for amplicon sequencing. While this marker has been used to assess some Actinobacteria groups from environmental samples, it has mostly been used with regard to particular groups, such as Nocardia (Vautrin et al., 2021) or Mycobacterium (van der Wielen et al., 2013). Therefore, considering these features, the hsp65 gene may be a good candidate for detailed assessments of the whole Actinobacteria community in extreme environments, such as caves. To make this possible, a reference database encompassing the Actinobacteria taxonomic group is required, which should be publicly available for taxonomic identification.
Since the 1960s, the Lascaux Cave has undergone several microbial outgrowth events until recently and has been subjected to (often unsuccessful) chemical treatments, which have resulted in the formation of various visual marks (including stains on walls) that have threatened the Paleolithic artwork (F. Bastian et al., 2010;Bontemps et al., 2021;Martin-Sanchez et al., 2015). Lascaux was closed to the public in 1963 to ensure cave conservation; despite this, microbial imbalance persisted and surface alterations kept forming, leading to visual marks on wall surfaces (Bontemps et al., 2021). Although the few microbiological studies conducted in this cave have suggested that Actinobacteria is a prevalent bacterial group (Alonso et al., 2018(Alonso et al., , 2019Martin-Sanchez et al., 2012), a more complete investigation of the Actinobacteria community is lacking, especially at the species level. Actinobacteria possess broad metabolic variability, and key functional traits may be related to precise species within the Actinobacteria; therefore, we hypothesized that such an assessment would better discriminate Actinobacteria communities according to ecological conditions. The objective of this study was to assess the usefulness of hsp65 to complement 16S rRNA gene sequence data using a metabarcoding approach for a deeper, species-level appraisal of the Actinobacteria community in caves, and to assess variation according to cave condition and human disturbance (Mammola, 2019). We compared communities belonging to the phylum Actinobacteria between pristine limestone caves (Reille and Mouflon Caves) and anthropized caves (Lascaux and Rouffignac Caves), and within Lascaux Cave between marked (wall visual marks) and unmarked areas in different rooms (i.e., Sas-1, Passage, Apse, and Diaclase). The second objective was to create an hsp65-based Actinobacteria database that is accessible to the scientific community.   Bastian et al., 2010;De La Rosa et al., 2017;Martin-Sanchez et al., 2012). Mouflon and Reille Caves remain in a pristine state (Alonso et al., 2019). The entrance of Mouflon Cave was locked shortly after discovery, whereas Reille is occasionally explored by seasoned speleologists.

| Study site description and sample collection
In Lascaux Cave, samples were taken from the walls of four distinct rooms representing different cave environments that are located progressively farther from the entrance in the following order: Sas-1, Passage, Apse, and Diaclase ( Figure 1b) Within each room, samples were collected from visual marks (when occurring) and unstained surfaces, as previously described (Alonso et al., 2018). Approximately 50 mg of wall material samples (3-6 replicates) were collected using sterile scalpels and placed in liquid nitrogen for transportation to the laboratory, where they were stored at −80°C until DNA extraction. Six samples from Rouffignac, six from Mouflon, five from Reille, and forty-one from Lascaux (five from Sas-1, six from Passage B, six from Passage IP, eighteen from Apse, and six from Diaclase) were used for subsequent metabarcoding analyses (Table A1).
Sampling was performed in accordance with the cave rules and regulations.
2.2 | DNA extraction and amplicon sequencing DNA was extracted using the FastDNA SPIN Kit for Soil (MP Biomedicals) following the manufacturer's instructions, and adapted to the low sample size, as previously described (Alonso et al., 2018).
The elution step was achieved using two volumes of 50 μl elution buffer per sample.

| Processing and analysis of sequence data
Paired sequence reads obtained from sequencing were merged using Fast Length Adjustment of Short reads (FLASh; Magoč & Salzberg, 2011) with a maximum of 25% mismatches in the overlapping region. The sequences were then filtered and aligned using reference alignments of hsp65 gene sequences (this study) or 16S rRNA gene sequences from the Silva database V132 (Quast et al., 2012). Chimeric sequences were removed using the integrated Vsearch tool (Rognes et al., 2016) according to the MiSeq standard operating procedure (MiSeq SOP, February 2018; Kozich et al., 2013) in Mothur v1.39.5 software (Schloss et al., 2009). Sequence libraries were taxonomically assigned in Mothur using the hsp65 reference database for hsp65 sequences and the recreated SEED database subset of the Silva Small Subunit rRNA Database, release 132 (Yilmaz et al., 2014), adapted for use in Mothur (https://mothur.org/w/ images/a/a4/Sylva.seed_v132.tgz), as the reference database for the 16S rRNA gene sequences.
The sequences of plastids and mitochondria, and those not classified in the domain Bacteria (from the 16S rRNA gene sequence library), as well as sequences identified as homologs for hsp65 (groEL, in the database renamed as GROESL; see the hsp65 reference database design section; from the hsp65 sequence library), were discarded. The sequence library was clustered into operational taxonomic units (OTUs) using the Uparse pipeline in Usearch v10.0.240 software (OTUs at a 97% cutoff; Edgar, 2013) and Mothur (OTUs at a 99% cutoff). Two cutoffs were applied, as the 97% cutoff is used as a standard threshold for estimating bacterial diversity for the 16S rRNA gene. Therefore, it might underestimate total diversity, especially that of the hsp65 marker, which has higher variability than the 16S rRNA gene marker (see Figure 2). The tighter 99% cutoff was used for a more reliable estimation of Actinobacteria diversity at the lowest taxonomic ranks, particularly for the hsp65 marker (Chen et al., 2013). The OTU table was further processed using tools implemented in Mothur. For 16S rRNA genes, only sequences corresponding to Actinobacteria were used.
To compare Actinobacteria among the caves, 25 samples for hsp65 and 24 samples for 16S rRNA genes were used (Table A1).
Apse and Diaclase (using unmarked areas) were selected as representative samples for Lascaux Cave, except for the analysis of the core microbiome (i.e., common OTUs present in all caves), for which all samples from Lascaux Cave samples were included. To examine differences among locations in Lascaux Cave, 33 samples for hsp65 (16 from marked areas and 17 from unmarked areas) and 38 samples for the 16S rRNA gene (20 from marked areas and 18 samples from unmarked areas) were used (Table A1). For the 16S rRNA gene marker, samples from the Sas-1 unmarked area were missing for the analysis, because only a low number of replicates were successfully sequenced. To obtain the highest number of replicates per sample, all samples were included in the analysis, even when they were not common to both markers. A total of 37 common samples for both markers were used to compare these markers, that is, rarefaction analysis, comparison of diversity indices, and several taxa recovered at different taxonomy ranks.
Rarefaction analysis was performed to analyze the richness of the Actinobacteria OTUs based on two cutoffs (i.e., 97% and 99%) as a function of sample number. The alpha diversity indices (richness, Chao-1; evenness, Simpson evenness; and diversity, inverse Simpson) were calculated in Mothur (97% and 99% OTU cutoff). Significant differences among samples were calculated by analysis of variance (ANOVA) with Tukey's post hoc tests (p < 0.05; for caves or locations in Lascaux Cave), Student's t-tests, and F tests (p < 0.05; for marked/unmarked areas) in the Past 4.02 software (Hammer et al., 2001). Pairwise comparisons were made with the Metastats method using Fisher's exact tests (White et al., 2009), and the linear discriminant analysis (LDA) effect size (LEfSe; Segata et al., 2011) was calculated to identify the Actinobacteria OTUs that differed significantly among the respective caves, locations, and marked/ unmarked areas in Lascaux Cave (99% OTU cutoffs). Metastats was also used to compare relative abundances of Streptomyces sequences among the marked/unmarked areas (hsp65 marker). Venn diagrams were created in Mothur to test the number of OTUs that encompassed the core microbiome for the different caves and rooms (i.e., Sas-1, Passage, Apse, and Diaclase) in Lascaux Cave (99% OTU cutoff).
Taxonomic classification of OTUs was performed with Mothur using the hsp65 and 16S rRNA gene SILVA V132 reference databases, and species-level identification with the hsp65 marker was verified by BLASTN in the online databases of the National Center for Biotechnology Information (NCBI; US National Library of Medicine). In parallel, a phylogenetic analysis based on hsp65 was performed using retrieved sequences from caves belonging to most representative OTUs and those from the GenBank database after using the nucleotide BLAST (BLASTN) tool. Retrieved sequences from GenBank were those with the highest similarity with selected OTUs. Multiple alignments were generated using ClustalX (Thompson et al., 2003). MEGA6 software was used to compute the maximum likelihood of neighbor-joining phylogenetic trees from the observed amino acid sequence divergence (%; Tamura et al., 2013). The taxonomical composition of the core microbiome was constructed using the Krona tool (Ondov et al., 2011). The number of taxa recovered with the common samples by hsp65 and 16S rRNA gene analyses at different taxonomic ranks was calculated and displayed using Venn diagrams.
Co-occurrence networks were used to predict how marked/  (Fruchterman & Reingold, 1991). The average degree, eigenvector centrality, and modularity were computed, and the minimum number of links was filtered to five (using a degree range filter).

| In silico analysis
Differences in the variability of the hsp65 and 16S rRNA gene partial sequences between pairs of Actinobacteria strains at different taxonomical levels (e.g., intraclass, intraorder, intrafamily, and intragenus) were determined using genomes retrieved from the Integrated Microbial Genomes and Microbiomes database (IMG, University of California, https://img.jgi.doe.gov/).
For the analysis, sequences from the available genomes of 47 species from the order Streptomycetales (Streptomyces, Kitakatospora, Table A2), and 58 species from the order Corynebacteriales (family Corynebacteriaceae: Corynebacterium; family Gordoniaceae: Gordonia; family Mycobacteriaceae: Mycobacterium; family Nocardiaceae: Nocardia, Rhodococcus; Table A2) were used. Sequences were aligned and trimmed to the amplified length of both gene markers in BioEdit 7.2.5 software (Hall et al., 2011). Pairwise distances between sequences were calculated in Mothur v. 1.39.5 software (Schloss et al., 2009).
Plots were created and the coefficient of determination (R 2 ) was calculated in Microsoft Excel (2016).
F I G U R E 2 (a) Rarefaction curves for 37 common samples of the hsp65 and 16S rRNA genes at the 97% and 99% OTU cutoffs. The X-axis denotes the number of samples, and the Y-axis denotes the number of OTUs. (b) Pairwise molecular distances between sequences of the hsp65 (Y-axes) and 16S rRNA gene (X-axes) genes among Actinobacteria species from different taxonomical levels with the correlation coefficient (R 2 ) for each equation. OTU, operational taxonomic unit 2.5 | hsp65 reference database design For taxonomic assignment using the hsp65 marker, the reference database of the Actinobacteria hsp65 gene was constructed similarly to that of the rpoB marker by (Ogier et al., 2019). Here, we built a new reference library based on the hsp65 gene encompassing both newly obtained and public sequences from the most important Actinobacteria taxa. Part of this new reference library was constructed with Nocardia and Gordonia sequences from all type strains, from strains from clinical collections (supplied by the French Observatory for Nocardiosis, OFN http://ofn.univ-lyon1.fr), and from environmental strains (supplied by CRB-EML http://eml-brc. org/) sequenced by our means for which sequence quality was verified by double sense sequencing. For the Mycobacterium genus, hsp65 sequences were available via the BIBI database (https:// umr5558-bibiserv.univ-lyon1.fr/lebibi/lebibi.cgi), which is fairly robust as it contains all described Mycobacterium species. The rest of the sequences was collected from the online IMG and NCBI databases using sequences homologous to the 65 kDa heat-shock protein (hsp65) of the reference strain Nocardia asteroides ATCC 14759 (Rodríguez-Nava et al., 2006) amplified by the TB11 and TB12 primers (Telenti et al., 1993). Paralogs of heat-shock proteins, such as groEL genes from Escherichia coli (Colaco & MacDougall, 2014;Duchêne et al., 1994; C. M. S. Kumar et al., 2015), were identified using maximum likelihood, FastTree 2.1 (Price et al., 2010), and were retained in the reference database (renamed GROESL sequences) as an outgroup that enabled amplified sequences belonging to this variant to be discarded. This database was named "ACTIhsp65_V1.0.0.fas" and is available at https://doi. org/10.5281/zenodo.5576073 The hsp65 database contained 198 genera and 1066 unique taxa, the whole represented by 5165 sequences, of which 1782 were denoted as GROESL at the time of the analysis. TB11 and TB12 primer complementarity with sequences from the reference database was evaluated using the OligoAnalyzer 3.1 tool (http:// www.idtdna.com/calc/analyzer). Percent coverage of primers for Actinobacteria genomes from different classes and orders was determined with primer BLAST using the RefSeq reference genomes (NCBI). Only targets with no mismatches for the last three nucleotides at the 3′ end and no more than four overall mismatches were included.

| RESULTS
3.1 | Sequence polymorphism of hsp65 versus 16S rRNA genes for cave Actinobacteria The rarefaction curves of 16S rRNA genes for both 97% and 99% cutoffs were more tilted, indicating an accumulation of identical OTUs due to their repeated sampling, which did not increase with a stricter cutoff, in contrast to hsp65 (Figure 2a). These results are consistent with the in silico analysis of hsp65 and 16S rRNA gene partial sequences obtained from 105 genomes (47 species from the order Streptomycetales and 58 species from the order Corynebacteriales) available on IMG. Pairwise distances between hsp65 sequences, even for closely related species, were higher than those between 16S rRNA genes ( Figure 2b).
The hsp65 primers targeted hsp65 sequences from the reference database with lower numbers of mismatches compared with the outgroup (i.e., distant homologs renamed in the database as GROESL, GroEL of E. coli, and Chloroflexi members; Table A3).
However, the primers only targeted select groups from the class Actinobacteria, as demonstrated in the analysis with representative genomes from online databases of the NCBI (US National Library of Medicine; Table A4). The number of taxa identified at class or lower taxonomic levels is shown in Figure A1. Data indicated that 20% of the classes, 32% of the orders, 36% of the families, and 25% of the genus-level taxa detected with hsp65 were also recovered through 16S rRNA gene analysis. Although there was a trend for a higher number of taxa uniquely detected by the hsp65 marker at lower taxonomic ranks, contrary to the 16S rRNA gene marker, hsp65 sequences detected a lower number of taxa overall. However, at the species level, s 168 species were recovered with the hsp65 marker only, and the highest species numbers were found within the Streptomyces, Mycobacterium, and Nocardia genera (Table A5).
In summary, the 16S rRNA gene marker was more suited for broad taxonomic profiling of the Actinobacteria community. In comparison, the hsp65 marker could help distinguish between selected Actinobacteria species, which the 16S rRNA gene marker failed to achieve; therefore, it complemented the results obtained with the 16S rRNA gene marker (Table A5).

| Actinobacteria diversity in caves based on hsp65 versus 16S rRNA genes
In the diversity analysis based on 37 common samples of both markers, the hsp65 marker recorded higher richness (Chao-1 index) and evenness (Simpson evenness index) than the 16S rRNA gene marker when considering the 99% OTU cutoff (Table A6).
The diversity analysis based on all samples indicated that Rouffignac Cave exhibited the highest Actinobacteria richness (Chao-1 index, for the hsp65 marker; Figure A2) but lower evenness (Simpson evenness index; for both markers; Figure A2) and diversity (inverse Simpson index; for the 16S rRNA gene marker; Figure A2) among other caves. In contrast, the Lascaux and Rouffignac Caves had lower evenness than the pristine Reille (for the 16S rRNA gene marker) and Mouflon (for the hsp65 marker; Figure A2) caves. When comparing locations within Lascaux Cave, significantly higher Actinobacteria richness was detected in Sas-1 compared with Passage B using only the 16S rRNA gene marker ( Figure A2). The average diversity of the unmarked areas was higher than that of the marked areas for both markers ( Figure A2).
3.3 | Actinobacteria community structure in caves based on hsp65 versus 16S rRNA genes NMDS analysis was conducted to compare the Actinobacteria communities in caves (97% OTU cutoff). These communities did not differ significantly between pristine caves, Reille and Mouflon (AMOVA), based on both hsp65 and 16S rRNA genes. NMDS distinguished three groups of communities corresponding to Lascaux, Rouffignac, and the two pristine caves, Reille and Mouflon, which was further supported by AMOVA (Table A7 A with both markers, except that differences between Passage IP, Sas-1 (both markers), and Diaclase (hsp65) were not significant due to the low number of samples (Table A7 A). Overall, the effect of surface alterations (i.e., marked vs. unmarked areas) in Lascaux Cave was significant for both markers (Table A7 A), with higher Actinobacteria homogeneity in unmarked versus marked areas (Table A7 B). Passage B presented the highest community stability compared with Passage IP, Apse, and Diaclase (Table A7 B), which was reflected in the NMDS plot, where the marked and unmarked areas in Passage IP were the least separated ( Figure 3).
3.4 | Core and cave-specific Actinobacteria microbiomes based on hsp65 versus 16S rRNA genes Venn diagrams showed that the number of OTUs shared between different caves (Rouffignac, Lascaux, Mouflon, and Reille) was two  Clusters of correlated OTUs dominated the unmarked areas, which also differed significantly in the Apse and Diaclase based on F I G U R E 6 Co-occurrence networks of Actinobacteria OTUs: (a) hsp65 marker and (b) 16S rRNA gene marker; OTUs differed significantly between marked (black) and unmarked (white) areas in Lascaux Cave, and those that did not differ between these areas (gray) using Metastats (p < 0.05). The letters indicate OTUs that were specific for the respective Lascaux Cave locations (S, Sas-1; PB, Passage banks; PI, Passage inclined planes; A, Apse; and D, Diaclase) using LEfSe (p < 0.03). Strong significant connections (Spearman's correlation >0.8 and <0.5 for 16S rRNA gene and <0.35 for hsp65) are displayed (99% OTU cutoff). Red lines indicate negative correlations and black lines indicate positive correlations. OTU, operational taxonomic unit both gene markers ( Figure 6). Additionally, OTU clusters that dominated the marked areas were typical for the Diaclase and Sas-1 areas based on the 16S rRNA gene marker. There was a correlation between OTUs from such distinct areas (Diaclase and Sas-1), suggesting that the marked-area factors affected the Actinobacteria communities in both locations ( Figure 6). However, the hsp65 gene marker revealed that the OTU clusters typical of marked areas were not typical of any Lascaux Cave locations. Finally, based on the 16S rRNA gene marker, the highest number of correlations between OTUs that did not differ between the marked and unmarked areas was found for those that dominated in Passage B, Sas-1, and Diaclase ( Figure 6). For those OTUs, the location factor was more influential than the marked/unmarked areas. Networks based on both markers

| DISCUSSION
This study represents an attempt to exploit the potential of hsp65 as a taxonomic marker to assess the whole Actinobacteria community in an extreme environment, such as caves, highlighting its advantages over the 16S rRNA gene marker. When comparing both genes, Venn results revealed higher variability for hsp65 sequences compared with Actinobacteria 16S rRNA gene sequences, and higher OTU richness (at the 99% cutoff) for Actinobacteria. Moreover, using the hsp65 marker, we found 168 species, which could not be identified in the same samples using the 16S rRNA gene marker. Therefore, the low interspecific polymorphism of 16S rRNA genes leads to an underestimation of diversity, which might be compensated for by using the hsp65 marker, as found for other protein-coding genes, such as rpoB (Vos et al., 2012). However, a higher number of higherrank taxa were identified with the 16S rRNA gene compared with the hsp65 marker, which highlights gaps in the hsp65 reference database, preventing the complete taxonomic assignment of sequences.
Another reason might be related to primer mismatches. Indeed, in silico analysis revealed that the primers amplified the gene from diverse Actinobacteria genomes; however, targets without perfect homology might be amplified with lower efficiency, leading to the underestimation of some taxa (Deagle et al., 2014). Overall, the data indicated that hsp65 is a suitable complement for the 16S rRNA gene marker in high-throughput sequencing methods for Actinobacteria, especially for taxonomic analysis at the species level; however, primer bias may need to be considered for community analyses at higher taxonomic ranks. Moreover, the hsp65 reference database is biased against representation in nature, due to uneven selection of restricted taxa, especially compared with those available in online databases. The 16S rRNA gene is available to the scientific community for all species because it must be mandatorily supplied as part of the submission. However, hsp65 gene submission is not necessary when describing a new species, explaining the lack of databases. The availability of bacterial genomes is increasing and in the future, many more sequences will be added that will increase the coverage of our database, allowing for more complete detection of taxa within the Actinobacteria.
In this study, the combination of hsp65 and 16S rRNA gene markers revealed that Actinobacteria diversity could vary according to anthropization in caves, on both higher (i.e., between caves) and lower (i.e., between locations within one cave) spatial scales. Indeed, the actinobacterial communities of anthropized (Lascaux and Rouffignac) and pristine caves (Reille and Mouflon) differed. The occurrence of these caves in the same region and limestone vein probably facilitated this observation, as geographic distance (Barraclough et al., 2012) and geological type (Zhu et al., 2019) are significant factors influencing cave communities. The number of Actinobacteria that could be identified as part of the core microbiome was relatively low, which suggested that different cave features can lead to different microbial communities (Alonso et al., 2019). Nevertheless, taxa belonging to this core microbiome (Pseudonocardiales, Streptomycetales, and Gaiellales) can resist anthropogenic disturbance (Shade et al., 2012).
Mouflon and Reille are pristine caves with evenly distributed taxa in the Actinobacteria community, possibly as a result of their stable environments without man-made disturbances (Mammola, 2019).
Typical pristine cave taxa that were identified included Gaiellales (Zhu et al., 2019) and Actinobacteria related to the environmental clone MB-A2-108_ge (De, 2019;, according to the 16S rRNA gene marker. High proportions of Streptomyces and unclassified Actinobacteria based on hsp65 data suggest that each marker identifies particular Actinobacteria community members. Less anthropized areas, such as the pristine caves of Reille and Mouflon, and the Diaclase area of Lascaux, yielded a higher proportion of "Actinobacteria unclassified" when using the hsp65 gene compared with the 16S rRNA gene in the same zones. Thus, these areas host species of families whose hsp65 gene sequences are missing in our database. Species may already be described for which the hsp65 sequence is not available, or they may simply represent new and as yet undescribed species. Therefore, our database, which was conceived to reach the species level, seems to be more suitable to identify sequences derived from anthropized environments. Environmental selection in oligotrophic habitats may promote rapid molecular evolution (Kuo & Ochman, 2009;Sagova-Mareckova et al., 2015) and result in a high number of novel species isolated from caves (De, 2019). A similar result was also found for the Streptomyces genus (Hamm et al., 2017), whose genomes may undergo a high rate of evolution (Cheng et al., 2015). For pristine caves, high proportions of Streptomyces were obtained from the Diaclase, which is an isolated compartment that has been exposed to little human impact within the anthropized Lascaux Cave. This observation seemed to confirm the hypothesis of Rangseekaew et al., who reported that, not only individual caves but also less-exposed locations within anthropized caves maintained typical Actinobacteria cave communities, which probably include novel taxa (Rangseekaew & Pathom-Aree, 2019).
Contrary to a previous study based on the whole bacterial community (Alonso et al., 2019), high richness but low evenness was observed in the visited Rouffignac Cave compared with the pristine caves, based on hsp65 analysis. The high richness of Actinobacteria in this cave might reflect a disturbance effect, which promoted the cohabitation of ecologically different microorganisms (Galand et al., 2016).
Other opportunistic pathogenic species such as N. farcinica were found in SAS-1, which is one of the most anthropized areas of Lascaux Cave. Its presence may be associated with human disturbance. The OTUs associated with this species were very close to the GenBank sequences obtained from a patient with pulmonary coinfection with M. tuberculosis (KF432743.1) and Madura foot (CP031418.1; Y. Y. Zhang et al., 2013). Other opportunistic pathogenic species of Streptomyces genera such as S. albus could also be found, at a low proportion, in Rouffignac cave, which may also be a consequence of human presence (Martín et al., 2004).
However, pathogenic species were also found in pristine caves, including N. cyriacigeorgica, which is one of the most frequent Nocardia opportunistic pathogenic species in France. The finding in Mouflon represents the first detection of this species in caves. Thus, this species may belong to the autochthonous core microbiome of the caves. The OTUs associated with this species were close to the GenBank sequences obtained from patients (EF127505.1, EF127507.1, and EF127504.1) and urban sediments (VBUT00000000.1) positioned in phylogroups I and III of the N. cyriacigeorgica complex (Schlaberg et al., 2008;Vautrin et al., 2021).
Our study has succeeded in positioning retrieved sequences in different phylogroups inside this complex.
As shown by Zoropogui et al. (2013), the whole genome of N.
cyriacigeorgica presents the acquisition of new genetic elements allowing it to thrive in multiple environments. This is confirmed by its isolation in a cave (this study), in the urban sediments of an infiltration basin (Vautrin et al., 2021), or even in humans. Its presence under different extreme environments could explain the evolution of this species into several different phylogroups, which may possess different levels of virulence. The sequences found in this cave may belong to the less virulent and most preserved phylogroups, unlike phylogroup II (the most virulent of this species), which was not identified here.
Other opportunistic pathogens belonging to the Mycobacteriaceae family were also identified in Lascaux Cave. The two most proportionally abundant species, M. lentiflavum and M. algericum, are causal agents of pulmonary infections in humans (Chida et al., 2021) and animals (Sahraoui et al., 2011), respectively. In addition, the presence of these species has been reported in aquatic mediums (Makovcova et al., 2015;Moradi et al., 2019;Mrlik et al., 2012).
OTUs obtained in our study that belonged to these two species were close to the GenBank sequences; one of them was obtained from Pinna nobilis (GenBank submission code MN854410.1). Thus, the presence of these species could be due to human input or the presence of water, possibly due to infiltration.
This study reports the occurrence of the most frequent opportunistic pathogenic Actinobacteria in France from caves. It is unknown whether environmental strains from these species present any pathogenic potential; however, the occurrence of these species in caves suggests a possible allochthonous input by tourists. In turn, this raises potential health concerns, as previously proposed , although Actinobacteria disease related to cave visitation has never been reported. In contrast, Jiangellaceae are a family typical of pristine caves (Rangseekaew & Pathom-Aree, 2019) and are another important group in Rouffignac Cave. According to the intermediate disturbance hypothesis (Roxburgh et al., 2004), the coexistence of taxa typical of pristine or anthropized conditions in the same habitat suggests that Rouffignac Cave is under intermediate The actinobacterial communities in caves were highly specific even within very short distances (Zhu et al., 2019) since each Lascaux room hosted its own Actinobacteria population. Lascaux Cave was typified by high relative abundances of Mycobacterium. Similar to Nocardia in Rouffignac, this group can signify external contamination Modra et al., 2017), although it is commonly found in natural environments (Kopecky et al., 2011). Confirming the previous results from Alonso et al. (2018), which were based on whole bacterial community analysis, the different geological substrates of Lascaux's Passage had a stronger role than the occurrence of visual marks; however, the marks had formed in the past and were stable. The effect of the geological substrate could also be observed when comparing Diaclase and Apse, because these two locations displayed different Actinobacteria communities, especially in unmarked areas, as shown by the network analysis. Relative isolation and distance from the cave entrance may represent efficient filters for alien microorganisms and enable the maintenance of the cave oligotrophic community, typical for unmarked areas Mammola, 2019).
Our results showed that the Actinobacteria community, mainly Streptomyces taxa, varied according to the presence of visual marks.
Marked areas, where Streptomycetaceae dominated, exhibited much lower diversity compared with unmarked zones, where Pseudonocardiaceae were prevalent. Unfortunately, comparisons of retrieved sequences with those in public databases yield low similarity scores that do not permit species associations to be made; therefore, we cannot hypothesize whether the presence of these OTUs is due to human disturbance. Similar to our study, the group Pseudonocardiaceae has been confirmed as a true cave rock dweller (Zhu et al., 2019). Moreover, this group has been identified in pigmented zones of pristine caves (Lavoie et al., 2017). Therefore, retrieved sequences may correspond to autochthonous species of the caves. In addition, members of the group Pseudonocardiaceae degrade complex molecules (Lavoie et al., 2017;Porca et al., 2012), which suggests that the unmarked areas were not completely oligotrophic (Tomczyk-Żak & Zielenkiewicz, 2016), perhaps as a result of past treatment with biocides. However, the higher community diversity in unmarked areas compared with marked areas suggests that cooperative relationships prevail, which is typical for oligotrophic cave environments (Tomczyk-Żak & Zielenkiewicz, 2016). In this study, specific features of visual marks were not considered according to room location and geological substrate, as the number of replicates was not sufficient; therefore, we focused on findings obtained at the scale of all visual marks considered together.
The lower diversity of marked areas might result from the competitive advantage of invading species (Hamm et al., 2017;Van Elsas et al., 2012). In isolated cave systems, invasions are mainly expected from the cave entrance, similar to Altamira Cave, where Streptomyces is the most dominant group, especially in stains (Groth et al., 1999). Our study reports the detection of S. mirabilis-like grouping, S. fulvissimus, and S. niveus, in caves. These species, which were found in stained areas of Lascaux Cave, have not been previously reported as pathogenic, but rather as producers of secondary metabolites with antimicrobial activity, some of which are used for bioremediation (El-Sayed, 2012;Flinspach et al., 2014;Kominek, 1972;Myronovskyi et al., 2013;Saberi-Riseh & Moradi-Pour, 2021;Schütze et al., 2014). Regarding S. fulvissimus and S. niveus, the OTUs belonging to these species were close to the GenBank sequences obtained from torrent muddy soil (CP071044.1), rhizosphere (CP048397.1), and deep-sea sediment (CP018047.1).
This seemed to confirm the ability of these species to survive in extreme environments away from human activity. Regarding the S. mirabilis-like grouping, the most proportionally abundant species in Lascaux stains found in this study, OTUs were close to the GenBank sequences obtained from fresh creek bank soil (CP074102.1). The set of OTUs associated with this group is positioned in a phylogenetic clade with only 60% bootstrap, which suggests that they could represent a new taxon within Actinobacteria communities, producers of stains in Lascaux Cave. Considering the reported antimicrobial activity of S. mirabilis, we believe that the ability of this probable new taxon to produce secondary metabolites may assist, through interactions with other microorganisms, in its survival in this stained environment.
Hypothetically, Streptomyces might also contribute to pigment production (Abdel-Haliem et al., 2013;Cuezva et al., 2012), or affect the pigment-forming fungi in Lascaux Cave (De La Rosa et al., 2017;Frey-Klett et al., 2011) based on their ability to cooperate with fungi and promote mycelial extension and secondary metabolite production (Frey-Klett et al., 2011). A previous study revealed the cooccurrence of Streptomyces and pigment-forming fungi, Acremonium and Exophiala, in Lascaux Cave, although only in unmarked areas (Alonso et al., 2018). Moreover, S. mirabilis is resistant to heavy metals (Schütze et al., 2014); this might play a role in the formation of pigmented marks since protection against toxicity of heavy metals and other chemicals has also been linked to melanin production (Nosanchuk & Casadevall, 2003). Therefore, the role of Streptomyces in the development of visual marks should be studied in more detail to evaluate whether this group may represent a missing link in our understanding of stain development in Lascaux Cave.

| CONCLUSION
Despite the limitations of our database, which yielded a high number of unclassified sequences, the hsp65 gene demonstrated its utility to complement the 16S rRNA gene which cannot resolve species but encompasses a larger part of the actinobacterial community. The use of hsp65 gene allowed us to detect, with a robust bootstrap, several species not previously reported in caves, demonstrating that it can be used reliably to differentiate Actinobacteria species in such environments. We expect these limitations to be addressed with the enrollment of new Actinobacteria sequences. This can be performed periodically by requesting other public databases, such as BLAST, or by enrolling sequences of previously identified strains isolated from environmental studies. The Actinobacteria community at different spatial levels reflected the natural quality of the caves and different locations within Lascaux Cave. Anthropization was shown to shape Actinobacteria communities, which altered typical cave communities and was exhibited by the Actinobacteria indicator taxa. Even if rare opportunistic pathogens from the Actinobacteria phylum were detected in anthropized caves, it is unknown whether these detected microorganisms represent an important health risk for visitors. We also found that marked areas on the surfaces of Lascaux Cave are important for shaping Actinobacteria communities, with particular reference to Streptomyces species, whose role in the formation of visual marks requires further research. Further studies are needed to isolate the Actinobacteria taxa detected in these caves to determine the true infectious risk associated with these pathogens and their role in the pigmentation of the marked areas. A better understanding of Actinobacteria functioning in caves will be important to guide efforts in the conservation of Paleolithic cave heritage.

CONFLICTS OF INTEREST
None declared. F I G U R E A5 Phylogenetic tree showing the taxonomic diversity and allocation of hsp65 sequences from caves retrieved by a metabarcoding approach. OTUs were used to build this figure. The evolutionary history was inferred using the neighbor-joining method (Saitou & Nei, 1987) and reveals the relationship of sequences from caves obtained in this study with sequences belonging to type and reference strains classified within the Actinobacteria class. The tree is based on the Kimura two-parameter method (Kimura, 1980) with the confidence values of the branches determined by bootstrap analyses (Felsenstein, 1985)  T A B L E A3 Percentage complementarity between the hsp65specific forward (TB11) and reverse (TB12) primers and the hsp65 sequence database (hsp65 sequences and GROESL sequences, paralogs to hsp65) with respect to the number of mismatches